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ABSTRACT 



Context. Radiative hydrodynamic simulations of solar and stellar surface convection have become an important tool for 
exploring the structure and gas dynamics in the envelopes and atmospheres of late-type stars and for improving our 
understanding of the formation of stellar spectra. 

Aims. We quantitatively compare results from three-dimensional, radiative hydrodynamic simulations of convection 
near the solar surface generated with three numerical codes (CD^BOLD, MURaM, and Stagger) and different simulation 
setups in order to investigate the level of similarity and to cross-validate the simulations. 

Methods. For all three simulations, we considered the average stratifications of various quantities (temperature, pres- 
sure, flow velocity, etc.) on surfaces of constant geometrical or optical depth, as well as their temporal and spatial 
fluctuations. We also compared observables, such as the spatially resolved patterns of the emerging intensity and of the 
vertical velocity at the solar optical surface as well as the center-to-limb variation of the continuum intensity at various 
wavelengths. 

Results. The depth proflles of the thermodynamical quantities and of the convective velocities as well as their spatial 
fluctuations agree quite well. Slight deviations can be understood in terms of differences in box size, spatial resolution 
and in the treatment of non-gray radiative transfer between the simulations. 

Conclusions. The results give confidence in the reliability of the results from comprehensive radiative hydrodynamic 
simulations. 

Key words. Methods: numerical - Sun: photosphere - convection - hydrodynamics - radiative transfer 



1. Introduction 



Comprehensive (magneto)hydrodynamic simulations have 
become an essential tool for studying near-surface con- 
vection in the Sun and other cool stars, together with 
the structure an d gas dynamics in their atmospheres (e.g., 
iNordlund et al.l . 2009). These simulations attempt to in- 
clude all relevant physics, such as three-dimensional (3D), 
time-dependent, compressible hydrodynamics, partial ion- 
ization and molecule formation as well as non-gray and non- 
local radiative transfer, in order to provide a 'realistic' rep- 
resentation of the physical stratification and macroscopic 
gas flows in the external stellar layers. For a direct com- 
parison with observations, spectral line profiles, continuum 
intensity and polarization maps are calculated on the ba- 
sis of the simulation results. This comparison serves as a 
means of validation of the simulations and also as a tool 
for interpreting the ob servational results in terms of basic 
physical quantities (cf. lUitenbroek fc Criscuoli l201lt) . 



Although various codes are now being used to per- 
form comprehensive simulations of solar and stellar (mag- 
neto) convection and an extensive body of simulation re- 
sults has already been published, so far no systematic at- 
tempt has been made to cross-validate codes by quantita- 
tively comparing numerical results. In this paper, we at- 
tempt to fill this gap, at least partially, and compare the 
solar models computed with CD^BDLD, MURaM, and Stagger, 
three independent and widely used 3D, radiative (mag- 
neto)hydrodynamic simulation codes. Apart from these 
codes, a number of other codes for 3D simulations of solar 
and stellar surface convection including (full or simplified) 
radiative transfer have been developed and utilized by var- 
lous groups (e.g..lStein fc Nordlundl. Il998t iRobinson et all 
2003; Heine mann et all. 120061: lAbbettl. 120071: lUstvugov . 
i2009t iMuthsam et all . 120101 : iGudiksen et al.l . l20lil l 



The purpose of our study is not a comparison of the 
numerical approaches of CO'^BOLD, MURaM, and Stagger per 
se, which would require using an identical setup in terms 
of box size, spatial resolution, and input material quanti- 
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Table 1. Numerical methods'^ used in the codes. 



Code 


spatial scheme temporal scheme 


RT scheme 


# rays 


# bins 


CO^BOLD 


Roe-type Riemann 


long characteristics 


17 


12 


MURaM 


4*-order FD 4"'-order RK 


short characteristics 


12 


4 


Stagger 


6*-order FD S'^'-order RK 


long characteristics 


9 


12 



^ FD: finite differences, RK: Runge-Kutta, RT: radiative transfer. 



ties such as opacities and equation of state. We rather wish 
to investigate how far simulations made for different appli- 
cations are consistent in the basic properties of the sim- 
ulated stellar atmosphere and uppermost convection zone. 
Examples of such properties are the average profiles of vari- 
ous quantities as a function of geometrical or optical depth. 
According to this rationale, we chose for comparison 'stan- 
dard' simulations of the near-surface layers of the Sun that 
were carried out by the participating groups for different 
purposes. The CO^BOLD and Stagger simulations provide a 
solar reference atmosphere (as part of a large grid of stellar 
models) for spectrum-synthesis calculations and abundance 
studies; they thus focus upon a good representation of the 
energy exchange by non-gray radiative transfer. The MURaM 
simulation, on the other hand, represents a non-magnetic 
comparison model for a set of magnetoconvection simula- 
tions to study fine-scale magnetic phenomena, for which 
high spatial resolution is crucial. The question we address 
here is: how much do the basic properties of the 'solar mod- 
els' resulting from simulations with different codes, differ- 
ent input quantities, different setup in terms of box size and 
spatial resolution, and even different top boundary condi- 
tions deviate from each other? If the differences turn out 
to be marginal, this result could then be taken as a kind of 
'cross-validation' of the codes and as a basis for confidence 
in the reliability of the simulations. 



2. Codes 

All three codes considered here treat the coupled 
time-dependent equations of compressible radiative 
(magneto)hydrodynamics and radiative transfer in a 
three-dimensional geometry and for a stratified, partially 
ionized medium. The energy exchange between radia- 
tion and matter is accounted for through solving the 
equation of radiative transfer under the assumption of 
local thermodynamic equilibrium (LTE) with a Planckian 
source function. To reduce the computational load, the 
wavelength-dependence of the radiative transfer is t reated 
with the method of o pacity binning dN ordlund, Il982t 
iLudwid . IT99I lSkarthenl . l200a iVogler et aI.L .2004) . A brief 
overview of the numerical methods used is given in Tabled] 
In all local 'box-in-the-star' setup is employed: the 

computational domain is a rectangular 3D box straddling 
in height the photosphere and the uppermost few Mm of 
the convection zoncQ The simulation boxes are sufficiently 
extended in the horizontal directions (6-9 Mm) to contain 
about 15-40 convection cells (granules) at any given time. 



thus providing a statistically useful sample of the near- 
surface layers of the Sun. Periodic boundary conditions 
are assumed in the horizontal directions (side boundaries) 
while the bottom boundary is open and allows free in- 
and outflow of fluid. A fixed entropy density is prescribed 
for the infiowing fluid at the lower boundary. It can be 
interpreted as the entropy of the deep, almost adiabatically 
stratified convective envelope and controls the effective 
temperature of the simulated atmosphere. In addition, the 
gas pressure is kept constant across the bottom boundary. 
The three codes differ somewhat in their treatment of the 
upper boundary conditions, as outlined more specifically 
below. 

Results from all three codes considered here already 
passed various 'reali ty checks' by comparison with obser- 



vation al data (e.g., Danilovic et al. . 20081: Pereira et al 
'2009a"l?; IWedemever-Bohm &: Rouppe van der Voort , 



2009; H irzberger et al.L l2010tl . 



^ However, by number of scale heights, the simulations cover 
about a third of the total pressure range of the convection zone. 



2.1. CO'' BOLD 

The CO^BOLD code uses a numerical scheme based on a 
finite-volume approach on a fixed Cartesian grid. Operator 
splitting separates the various (usually explicit) operators: 
the (magneto)hydrodynamics, the tensor viscosity, the radi- 
ation transport, and optional source steps. Directional split- 
ting reduces the multi-dimensional hydrodynamics prob- 
lem to a sequence of ID steps. The advection step is per- 
formed by an approximate Riemann solver of Roe type, 
modified to account for a realistic equation of state, a non- 
equidistant grid, and the presence of source terms due to an 
external gravity field. Optionally, a 3D tensor viscosity can 
be activated for improved stability in extreme situations. 
Parallelization of CD'''BOLD is achieved with OpenMP. 

The top boundary condition provides transmission of 
waves of arbitrary amplitude, including shocks: typically, 
two layers of ghost cells are introduced, where the velocity 
components and the internal energy are kept constant and 
the density decreases exponentially with a scale height set 
to a controllable fraction of the local hydrostatic pressure 
scale height. This gives the possibility to minimize the mean 
mass flux through the open top boundary. 

The radiative transfer in the CO^BDLD simulation consid- 
ered here was computed using a long-characteristics scheme 
for rays with four inclination angles and four azimuthal 
angles plus the vertical, i.e. 17 rays in total. The values 
of the inclinations (/i = cos6'=1.000, 0.920, 0.739, 0.478, 
0.165) correspond to the positive nodes of the 10*^-order 
Lobatto quadrature formula, the four azimuthal angles co- 
incide with the grid directions (0=0, 7r/2, tt, and 37r/2). 
For the opacity binning, tables were constructed from a 
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Table 2. Optical depth ranges and wavelength ranges'* 
the opacity bins. 

CO'^BOLD 



of 



Bin 


log Tl 


log T2 > 


1 [nm] 


A 2 [nm] 


I 


U.-LU 




n 

u. 


oou. 


2 


0.15 


99.00 


550. 


100000. 


3 


0.00 


0.15 


0. 


600. 


4 


0.00 


0.15 


600. 


100000. 


5 


-0.75 


0.00 


0. 


650. 


6 


-0.75 


0.00 


650. 


100000. 


7 


-1.50 


-0.75 


0. 


100000. 


8 


-2.25 


-1.50 


0. 


100000. 


9 


-3.00 


-2.25 


0. 


100000. 


10 


-3.75 


-3.00 


0. 


100000. 


11 


-4.50 


-3.75 


0. 


100000. 


12 


-99.00 


-4.50 


0. 


100000. 


MURciM 


1 


0.00 


99. 






2 


-2.00 


0.00 






3 


-4.00 


-2.00 






4 


-99. 


-4.00 






Stagger 


1 


-1.46 


9.00 


0. 


380.9 


2 


-3.81 


-1.46 


0. 


380.9 


3 


-15.00 


-3.81 


0. 


380.9 


4 


-0.62 


9.00 


380.9 


562.4 


5 


-0.62 


9.00 


562.4 


2161.2 


6 


-1.50 


-0.62 


380.9 


642.6 


7 


-2.28 


-1.50 


380.9 


710.9 


8 


-10.00 


-2.28 


380.9 


1646.5 


9 


-0.62 


9.00 


2161.2 


100000. 


10 


-1.50 


-0.62 


642.6 


100000. 


11 


-2.28 


-1.50 


710.9 


100000. 


12 


-10.00 


-2.28 


1646.5 


100000. 



The bins in the MURelM simulation are wavelength- 
independent. 



timization was not necessary because a tabulated EOS is 
used. 

The CO^BOLD simu l ations considered in this paper was 
used by ICaffau et al.l (|2008[ ) for the determination of the 
solar thorium and hafnium abundances, and for subsequent 
studies of CNO and other element s. More detail s abou t 
the CD^BOLD code can be found in iFrevtag et al.' ('20021), 



fWedemever et al.l (|2004[) , and iFrevtag et al. (2008,, .2011.) 



2.2. MURaM 



The MURaM code ("Vogleil l2003t IVoeler et all I2005D uses a 
4*'*-order central difference scheme in space and a 4*''-order 
Runge-Kutta scheme for time-stepping. Arti ficial diffusivi- 
ties ar e treated with the scheme described in lRemoel et al.l 
()2009[ ). An open-top boundary condition is also imple- 
mented in the MURaM code, but for the simulation consid- 
ered here a stress-free, closed top (zero vertical velocity) was 
chosen in order to study how far this affects the mean strat- 
ification. MURaM uses the Message Passing Interface (MPI) 
framework for parallelization. 

Radiative transfer in the MURaM code is calculated with 
the short-characteristics method (jKunasz fc Aued . 1198811 
with bilinear interpolation. The angu lar integration is car- 
ried out according to the A4 scheme of lCarlsonl ()1963l) along 
three directions per o ctant, which corre sponds to 12 com- 
plete rays in total (cf. lBruls et al.l . [T999D . The opacity bin- 
ning for the non-gray radiative transfer is based on the 
opacity distrib ution functions from the ATLAS9 package 
()Kuruczj . Il993f ) and uses 4 bins. The thresholds in optical 
depth (see Table [5]) for the binning procedure were chosen 
in terms of log r , which is a hybrid of the Rosseland mean 
in the deeper layers and the Planck mean in the u pper lay- 
ers, with a smooth transi tion centered at r = 0.35 ( jLudw ij. 
119921 : IVoeler et al.l . l200l ). The EOS tables used for the sim- 
ulation considered her e are b ased on tables from the OPAL 
project (jRogers et all 19961) for a so l ar ga s mixture with 
abundances from Anders fc Greves"s3 (jl989t ). 



data set of MARCS raw opacities provided by B. Plez 
(priv. comm.; see also iGustafsson et al.l . I2008D , comprising 
continuous and sampled atomic and molecular line opaci- 
ties as functions of temperature and gas pressure at more 
than 10^ waveleng th points. The a dopte d chemical compo- 
sition comes from lAsplund et al.l ([200^ . Each wavelength 
of the original opacity sampling data was sorted into one 
of twelve representative bins, according to wavelength and 
Rosseland optical depth where the monochromatic optical 
depth unity is reached in a ID standard solar atmosphere. 
The thresholds for the opacity bins used for the present 
CO^BDLD solar simulation are given in Table (H For each 
opacity bin, the tabulated opacity is a hybrid of Rosseland 
and Planck means over all frequencies of the bin, such that 
it approaches the Rosseland mean at high values of the 
optical depth, and the Planck mean at low values, with a 
smooth transition centered at Rosseland optical depth 0.35. 
Assuming LTE and pure absorption, the source function in 
each bin is computed as the Planck function integrated over 
the frequencies associated with the respective bin. 

T he equ ation of state (EOS) used in CO^BOLD follows 

IWolj (jl983f) . It accounts for the partial ionization of hy- 
drogen and helium, as well as for H2 molecule formation. 
In contrast to Wolf's approach, all pressure-temperature 
regions are treated homogeneously since performance op- 



2.3. Stagger 

The StaRger cod e, [originally developed by 

iGalsgaard & Nordlundl . [19930, uses a 6"*-order finite 
difference scheme in space with 5*^-order interpolations. 
Scalar variables (density, internal energy, and temperature) 
are volume-centered, while momenta are face-centered. 
The hydrodynamic variables are advanced forward in time 
using a 3'''^-order Runge-Kutta scheme. Boundaries are 
periodic horizontally and open vertically, b oth at the top 
and a t the bottom. The EOS is taken from iMihalas et all 
(1988^ and accounts for the effects of excitation, ionization, 
and dissociation of the 15 most abundant elements and of 
the H2 and H2 molecules. Parallelization of the Stagger 
code is carried out via MPI. 

The radia tive-trans f er eq uation is solved with a 
Feautrier-like ([Feautrieil . [1964) scheme along eight inclined 
rays (two inclination angles, four azimuth angles) plus the 
vertical, and using an opacity-binning scheme with twelve 
bins for the frequency dependence. The total radiative heat- 
ing rate at the center of each grid cell is computed by 
adding the partial contributions from each direction and 
opacity bin with the appropriate weight. The values of the 

^ see also http : //www . astro .ku.dk/ ~kg/Papers/MHD_code .ps .gz 
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inclination angles (/i=cos 6'=0.155, 0.645, 1.000) and their 
associated weights (w^=0.376, 0.512, 0.111) correspond to 
the nodes and weights of the 3'^'^-order Radau quadrature 
formula; the four azimuthal angles are equidistant (0=0, 
7r/2, TT, and 37r/2) and have equal weights. 



Table 3. Parameters of the simulation runs. 



Code 


Box [Mm^] 


h [Mraf 


grid resolution [km] 

{x,y,z) 


CO'^BOLD 


5.6 X 5.6 X 2.3 


0.88 


40 X 40 X 15.1 


MURaM 


9x9x3 


1 


17.6 X 17.6 X 10 


Stagger 


6 X 6 X 3.6 


0.88 


25.1 X 25.1 X 7.. .32 



X, y: coordinates in the horizontal directions; 
z: depth coordinate. 
^ height of top boundary above t — 1 



Table 4. Global properties of the simulated solar models. 



Code 




{5I/I)^oi[%] 


{5I/I)M%] 


CO'^BOLD 


5782.1 ± 12.6 


14.4 ±0.6 


21.8 ± 0.8 


MURaM 


5768.4 ±9.9 


15.4 ±0.3 


21.8 ± 0.3 


Stagger 


5778.4 ± 15.8 


15.1 ± 0.5 


22.1 ±0.8 



For the Stagger simulatio n considered here, the co ntin- 
uous opacity data came from iGustafsson et al.l (|1975[ ) and 
Trampedach (private communication), while sampled line 
opacities were taken from the M ARCS package (B. Plez, 
private communication; see also IGustafsson ct al., 200^. 
The adopted chemical comp osition for the s imula tion con- 
sidered here was taken from lAsplund et all ([2003). 

The opacity binning procedure implemented in the 
Stagger code is essentially based on the formulation by 
[Skartlici^ (j2000l) . Opacities are sorted into bins according to 
their wavelength and strength. As a measure of the opacity 
strength at a given wavelength, the Rosseland optical depth 
of formation of that particular wavelength is used. More 
precisely, the formation depth is defined as the point where 
the monochromatic optical depth in the vertical direction 
equals unity in a one-dimensional model constructed by 
taking the mean temperature-density stratification from a 
solar simulation. The thresholds in Rosseland optical depth 
and wavelength for the determination of bin membership 
are given in Table [21 Within each opacity bin, opacities 
are averaged and the source function contributions at the 
various wavelengths belonging to the bin are integrated. 
Mean-intensity-weighted average opacities and Rosseland- 
like mean opacities are adopted in the optically thin and op- 
tically thick layers, respectively. A bridging function is used 
for a smooth transition between the two averages near the 
optical surface. In the Stagger simulation considered here, 
the Planck function at the local temperature was chosen 
as the source function and the contribution of scattering to 
the total opacity in the optically thin layers was neglected. 



The simulation belongs to the s e ries th at was used in the re- 
cent analysis bvlAsplund et al] (|2009[ ) for the spectroscopic 
determination of solar abundances. A more comprehensive 
description of the opacity binning implementation and of 
the approximations i nvolved in current Stagger code sim- 
ulations are given in iCoUet et all (l20Tll) . 



3. Simulation runs and quantities for comparison 

The three codes were used to carry out simulations of 
near-surface solar convection without magnetic field. As ex- 
plained in the introduction, the simulation setups were dif- 
ferent, corresponding to the different research topics that 
the participating groups focus upon. Therefore, the nu- 
merical setups differ in terms of (horizontal and vertical) 
box size, grid resolution, number of opacity bins, and other 
features such as the top boundary condition. In all cases, 
the simulation boxes include the photosphere (up to about 
1 Mm above the optical surface) and the uppermost layers 
of the convection zone (between 1.4 Mm and 3 Mm below 
the optical surface, depending on the simulation). Table |3] 
gives various parameters of the simulation runs. Note that 
the Stagger code uses a non-equidistant grid of 230 cells 
in the vertical direction, with spacings ranging from 7 km 
around the optical surface and 32 km in the deepest parts 
of the simulation box. 

The simulations were run for several hours of solar time 
to reach a statistically stationary, thermally relaxed state. 
Nineteen snapshots taken at regular intervals and spanning 
in total a period of about two hours were considered for the 
analysis of each simulation. This choice was made to ensure 
that the effects of the 5-minute p-mode oscillations in the 
simulations are averaged out in the temporal means. 

The physical quantities considered for the compari- 
son are temperature, gas pressure, and turbulent pressure 
{pvl), as well as the vertical and horizontal velocity compo- 
nents. To obtain mean profiles as functions of depth, z, for 
each of these quantities, q{xi,yj, Zk,ti) = qijk,i, and their 
squares at the grid cells (xi,yj, Zk) and at time t = U, the 
averages over horizontal planes {zk = const.) were deter- 
mined, viz. 

^k,i = rr-im^yM (1) 

= ^EE^k/ ' (2) 

^ 2—1 J — I 

where rix and ny are the number of grid cells in the horizon- 
tal directions. Similarly, averages over surfaces of constant 
optical depth were determined by first calculating the op- 
tical depth along vertical lines of sight, T500, for the contin- 
uum opacity at 500 nm wavelength and then considering the 
quantities at fixed levels in the range —4 < logrsoo < 4. We 
also considered averages based on optical depth correspond- 
ing to the Rosseland mean opacity and found the results to 
be not significantly different from those with 500 nm con- 
tinuum opacity, so that we restrict ourselves to the latter 
case. 
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Fig. 1. Vertically emerging continuum intensity at 500 nm for single snapshots from the CO^BOLD (left), Stagger (middle), 
and MURaM (right) runs, drawn to scale. The gray scales cover, from black to white, the ranges 0.59-1.53 (CO^BQLD), 0.50- 
1.53 (Stagger), and 0.49-1.66 (MURaJ^) of the intensity normalized to the respective horizontal average. Axis units are 
Mm. 




Fig. 2. Vertical velocity at the average geometrical depth level of the surface r5oo = 1 for the same snapshots as in Fig.[TJ 
Downflows are shown in red, upflows in blue. The color table covers the range ±7 km s^^ in all cases; speeds outside this 
range are saturated. Axis units are Mm. 



From the vertical profiles Qj. ; and g^^. ; we determine and the temporal average of the spatial root- mean-square 
temporal averages over the = 19 snapshots considered (RMS) fluctuation. 



for each simulation, 

N 



1 

N 



(3) 



1=1 



the standard deviation among the snapshots, 

1/2 

1 ^ . 



N 



N ^ 



(9RMs)fc 



1 ^ 



9 kA 



{ikjY 



1/2 



(5) 



The index k in Eqs. ([SHSl) refers either to the geometrical 
depth level, z^, for averages over planes of constant geomet- 
rical depth or to the optical depth level, T5oo,fc, for averages 
over surfaces of constant optical depth. Strictly speaking, 
the latter averages are taken over the projections of the cor- 
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rugated surfaces of constant optical depth on a horizontal 
plane. 

4. Results 

Table H] shows the effective temperatures of the various sim- 
ulation models together with the disk-center bolometric and 
monochromatic continuum intensity contrasts at 500 nm. 
The standard deviations indicate the variability among the 
19 snapshots from each simulation run used in the analy- 
sis. The effective temperatures differ by about 14 K at most 
and the intensity contrasts agree fairly well with each other. 
In comparison to the other simulations, the variations from 
snapshot to snapshot are smaller in the MURaM case. This is 
probably because of the larger horizontal extension of the 
computational box. 

Unless stated otherwise, all quantities discussed in the 
following subsections refer to averages over the 19 snapshots 
from each of the simulation runs. 

4.1. Surface maps and histograms 

Figure [T] shows maps of the vertically emerging (disk- 
center) continuum intensity at 500 nm for snapshots from 
the three simulation runs. Fig. [5] gives the corresponding 
maps of the vertical velocity at the average geometrical 
depth level of the surface T500 = 1, where T500 is the con- 
tinuum optical depth at 500 nm wavelength. The runs with 
higher spatial resolution naturally show more small-scale 
details, but the basic structure and average size of the gran- 
ules and the correlation between the brightness and velocity 
are very similar in all simulations. The visual impression is 
confirmed by the similarity of the histograms of intensity 
and vertical velocity given in Fig. |21 On a logarithmic scale 
(right panels), the difference in spatial resolution becomes 
apparent at the extreme values, but otherwise there are no 
significant differences between the distributions. 

4.2. Mean stratification 

The upper panels of Fig. U show the geometrical-depth 
profiles of the horizontally averaged temperature and their 
(temporal) standard deviations (see Eq. H]) , the latter indi- 
cating the level of fluctuations of the mean profiles among 
the 19 snapshots used from each of the three simulation 
runs. The depth scales of the three models were aligned such 
that z = always refers to the average depth of the surface 
''500 = 1- The weaker fluctuations between the MURaM snap- 
shots compared to the other models are probably a result of 
the horizontally more extended computational box and the 
explicit damping of the fundamental box oscillation mode 
in this simulation. 

The absolute and relative differences between the mean 
temperature profiles from the three simulations are given 
in the lower panels of Fig. S) The colored bands in this 
plot (and in similar subsequent figures) represent the sam- 
ple standard deviation^ of the differences between pairs of 
models on the basis of the 19 snapshots from each simula- 
tion, indicating the range of scatter caused by the temporal 
variability of the averages in the relatively small simulation 

^ The standard deviation of the differences is equal to the 
square root of the sum of the variances according to the 19 
snapshots of the individual models. 



xCO^BOLD ■ 




z[Mm] 

Fig. 6. Horizontally averaged ratio of turbulent pressure to 
gas pressure. 



boxes. The temperature profiles agree fairly well, with rel- 
ative differences below 2% and overlapping bands of the 
standard deviation over most of the height range in the 
lower left panel of Fig. 2] The positive and negative ex- 
cursions near z — result from slight differences in the 
depth location and slope of the steep temperature gradient 
near optical depth unity. These differences can easily arise 
given the strong temperature dependence of the continuum 
opacity around T(r5oo = 1) ~ 6400 K and the differing 
vertical resolution of the simulations. In the photosphere 
(-500 km< z < 0km), the CO^BOLD and Stagger profiles 
deviate by less than 20 K from each other while the MURaM 
temperatures differ somewhat more, up to 60 K (i.e., at a 
level of about 1%). 

Figure [5] shows the depth profiles of gas pressure, Pgas, 
and of turbulent pressure, {pvj,), together with the respec- 
tive relative differences between the simulations. The ratio 
of turbulent to gas pressure is given in Fig. [51 The turbu- 
lent pressure reaches nearly 20% of the gas pressure slightly 
below the optical surface (where the vertical velocity fluc- 
tuations peak) and again about this value in the top layers. 

The relative differences of the gas pressure between the 
models are somewhat larger than those of the tempera- 
ture, especially in the uppermost layers. While the roughly 
depth-independent deviations in the layers below z = can 
be simply explained by a constant relative shift of the re- 
spective geometrical height scales, the bigger deviations in 
the photosphere are caused by the significant differences 
of the turbulent pressure in these layers between the sim- 
ulations (cf . the lower right panel of Fig. [5]) . To illustrate 
the effect, consider the simple case of an isothermal atmo- 
sphere and constant turbulent speed, Vz- The scale height 
of the gas pressure is then given by Hp = {c^-\-v1)/g, where 
Cs is the sound speed and g is the (constant) gravitational 
acceleration. If the turbulent pressure differs between two 
stratifications, the effect is cumulative: for instance, a dif- 
ference of 5% over 6 scale heights adds up to 30% of a scale 
height, leading to a significant pressure deviation in the 
upper layers. The higher turbulent speeds of the Stagger 
model imply a larger scale height and, therefore, higher 
pressure and density in the upper layers. In terms of opti- 
cal depth (relevant for observations), the deviations of the 
pressure stratifications are significantly smaller (cf. Fig.[TT|l 
since higher values of pressure and density lead to an up- 
ward shift of the iso-r surfaces (and vice versa). 
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Fig. 3. Linear (left) and logarithmie (right) histograms of the vertically emerging continuum intensity at 500 nm (upper 
panels) and the vertical velocity on the average height level of the surface r5oo = 1 (lower panels) averaged over all 19 
snapshots from each simulation (black: CO^BOLD, red: MURaM, blue: Stagger). Positive velocities correspond to upflows. 
Thirty bins were used in all cases. Each histogram was normalized such that the sum of the density function over the 
bins becomes unity. 



Figure [7] shows the relative RMS fluctuations of temper- 
ature and pressure, respectively, on surfaces of constant ge- 
ometrical depth (see Eq. [S]) . The temperature fluctuations 
(left panel) show a sharp maximum near the optical surface, 
where radiative cooling leads to strong temperature differ- 
ences between granular upflows and intergranular downflow 
regions. The development of shocks in the uppermost layers 
of the simulation boxes results in a second peak of the tem- 
perature fluctuations. In contrast, the pressure fluctuations 
grow monotonically outward and reach their maximum at 
the top of the simulated regions. 

The RMS values of the velocity components are shown 
in Fig. [HI The RMS of the vertical velocity (left panel) peaks 
near the optical surface, owing to the braking of the up- 
flows and the acceleration of the cool downflows as a re- 
sult of radiative cooling. Because the scale height decreases 
rapidly, most of the rising fluid that reaches the surface has 
to overturn very near to optical depth unity, so that the 
RMS value of the horizontal velocity (right panel of Fig. [5]) 
peaks only slightly higher than those of the vertical velocity. 
The RMS of both velocity components grow again in the 
upper, shock-dominated top layers of the simulation boxes. 
The somewhat lower RMS values of the CC'^BOLD model 
are probably related to the shallower computational box in 
comparison to the other simulations. The drop near the bot- 



tom boundary in the case of the MUReiM simulation is caused 
by a narrow layer of enhanced viscosity, which was intro- 
duced on grounds of numerical stability. In spite of these 
differences, all three codes show excellent agreement in the 
observable photospheric layers (—0.5 Mm < z < OMni). 

4.3. Photospheric structure 

For all observables originating in the solar photosphere, the 
profiles of the physical quantities as functions of optical 
depth are relevant. Since the surfaces of constant optical 
depth are not flat but strongly corrugated in the photo- 
sphere owing to granulation, the profiles of quantities av- 
eraged over surfaces of constant optical depth generally do 
not simply correspond to stretched or shifted profiles of 
horizontally averaged quantities. 

Figure [9] illustrates the relation between the geometri- 
cal and the optical depth scales. The left panel shows the 
mean geometrical depths of the iso-r surfaces as a function 
of continuum optical depth at 500 nm. The profiles are very 
similar: the maximum differences between the models are 
on the order of the vertical distances of the grid cells. The 
(spatial) RMS fluctuations of the depth of the iso-r surfaces 
are shown in the left panel of Fig. [HI They quantify the 'cor- 
rugation' of the iso-T surfaces, which reaches a maximum 
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Fig. 4. Upper panels: horizontally averaged temperature (left) and standard deviation of the mean temperature profiles 
corresponding to the 19 simulation snapshots (right) as functions of geometrical depth (z = 0: average depth of T500 = !)■ 
Lower panels: absolute (left) and relative (right) mean temperature differences between the models. 
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Fig. 5. Upper panels: horizontally averaged gas pressure (left) and relative differences between the models (right) as 
functions of geometrical depth. Lower panels: same for the horizontally averaged turbulent pressure, pturb = pui- 
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Fig. 7. RMS fluctuations of temperature (left) and pressure (right) on surfaces of constant geometrical depth. 





Fig. 8. RMS of the vertical (left) and horizontal (right) velocity on surfaces of constant geometrical depth. 
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Fig. 9. Left: geometrical depth averaged over iso-r surfaces (left) as functions of continuum optical depth T500. Right: 
spatial RMS fluctuations of the geometrical depth of the iso-r surfaces. 



of ^80 km somewhat below the optical surface, presumably 
because of the dominant effect of the cool downflow regions 
and the strong (positive) temperature dependence of the 
H~ continuum opacity. 

The profiles of temperature averaged over iso-r surfaces 
and the differences between the models are shown in Fig.lTOl 
In the photosphere, the biggest difference between the mod- 
els appears in the range —1.5 < logrsoo < —0.5, where 



the MURaM model is up to 40 K cooler than the Stagger 
model and 65 K cooler (~ 1%) than the CO^BOLD model. 
This probably is a result of the less detailed opacity bin- 
ning procedure used in the MURaM simulation (only 4 opacity 
bins compared to 12 in the other simulations). The devia- 
tion leads to a somewhat lower temperature gradient in the 
MURaM case between —2.5 < logrsoo < —1.5. In a narrow 
layer around logrsoo — 0.7 below the optical surface, the 
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Fig. 10. Temperature averaged over iso-r surfaces (left) and temperature difference between the models (right) as a 
function of optical depth T500 in the photospheric layers. 
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Fig. 11. Gas pressure profiles (left) and relative differences (right) averaged over iso-r surfaces as functions of optical 
depth T500. 
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Fig. 12. Relative RMS fluctuations of temperature (left) and pressure (right) on surfaces of constant optical depth tsoo- 



MURaM model is up to 200 K hotter than the other two mod- models in the photosphere generally amount to a few per- 
els. This could be related to the higher spatial resolution of cent. The biggest differences of up to 6.5% arise between 



the MURaM simulation, which leads to less artificial diffusion 
of heat between the hot upflows and cool downflows, thus 
maintaining bigger temperature fluctuations (cf. Fig. [13 



The profiles of gas pressure averaged over iso-r surfaces 
are shown in Fig. [TT] The relative differences between the 



the Stagger and MURaM models near logrsoo — —1. As ex- 
plained in the preceding subsection, these differences are 
related to the effect of the turbulent pressure on the pres- 
sure scale height. 
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Figure [T2] shows the spatial RMS fluctuations of tem- 
perature and pressure on iso-r surfaces. In the range —2 < 
logTsoo < 0, the horizontahy averaged temperature fluc- 
tuations are similar for all models. In the higher layers, 
the models deviate more strongly from each other. This is 
a result of the differences in the velocity fluctuations and 
turbulent pressure discussed above. 

The RMS values of the vertical and horizontal compo- 
nents of the fluid velocity are shown in Fig. [T31 The RMS 
values of both components in the photosphere are very sim- 
ilar for all three simulations, while the deviations become 
somewhat larger in the layers above. 

4.4. Center-to-limb variation of continuum intensity 

In order to estimate how much the differences between 
the models affect observable quantities, we considered the 
center-to-limb variation (CLV) of the monochromatic con- 
tinuum intensity. The intensities were calculated at four 
wavelengths (400 nm, 600 nm, 800 nm, and 2000 nm) along 
rays with 10 inclination angles, with fi — cos 9 ranging from 
/i = 1 (disk center) to /i = 0.1 (limb). For each inclina- 
tion and wavelength, the intensities were spatially averaged 
over the simulated surface areas, temporarily averaged over 
simulation snapshots, and azimuthally averaged over four 
(CD^BOLD, MURaM) or 12 (Stagger) equidistant azimuthal 
directions. 

The CLVs for the CD^BOLD simulation were computed 
with the spectrum synthesis code LinforSlQ. The NLTE 
EOS used in LinforSD is more detailed than that employed 
in CO^BOLD, because it has to provide the electron pres- 
sure and the number density for all individual atoms and 
ions. Likewise, the continuum opacities used in LinforSD 
are not fully consistent with the raw opacities from which 
the binned opacities in the CD^BDLD simulations are con- 
structed. However, both o pacities are based on the same 
chemical abundance mix ijCrevesse fc Sauvall . I1998L with 
the exception of CNO, for which the values A(C) = 8.41, 
A(N) = 7.80, A(0) = 8.67 are adopted). The emergent 
continuum intensity is obtained by integrating the trans- 
fer equation on a grid that is refined with respect to the 
original hydrodynamics grid, ensuring that the resolution 
in vertical optical depth is about Atross ~ 0.1. 

The CLVs for the snapshots of the Stagger simula- 
tion were computed with the line formation code SCATE 
(|Havek et al.l . 1201 Ih . using the same Feautrier-like long- 
characteristics solver as in the original simulation. SCATE 
employs the same EOS as the Stagger code simulation and 
the same continuous opacities used for the opacity binning. 

For the MURaM simulation, the CLVs were calculated us- 
ing a long-characteristics scheme with automatic grid re- 
finement in case of steep gradients in optical depth along 
the ray. It uses the same EOS as the MURaM simulation and 
the same continuum opacities as employed for the opacity 
binning (based on the opa city distributio n functions from 
the ATLAS9 package, see lKurtic3 . Il993l) . The continuum 
opacities are averages over 20 nm (for the CLVs at 400 nm, 
600 nm, and 800 nm) or over 90 nm (for 2000 nm). 

Fig. HH shows the resulting profiles for the three simu- 
lations and the differences between them. The profiles are 
very similar, the agreement between CO^BDLD and Stagger 
simulations being somewhat better (except at 400 nm) 

* http;/ /www. aip.de/~mst/Linfor3D/linfor_3DjTianual.pdf 



than that between MURaM and the other simulations. The 
somewhat stronger limb darkening mainly results from the 
slightly steeper temperature gradient in the lower photo- 
sphere of the MURaM model (cf. Fig. ITUl). 

5. Conclusions 

Although the three numerical solar models considered here 
(CO^BDLD, MURaM, and Stagger) result from codes with dif- 
ferent numerical methods and from simulation boxes of dif- 
ferent size and spatial grid resolution, their overall agree- 
ment is very good and encouraging. This does not only 
concern the mean (optical and geometrical) depth profiles 
of the basic quantities, but also the spatial fluctuations of 
these quantities as well as histograms of velocity and inten- 
sity, i.e., the dynamics and spatial structure of the simu- 
lated atmospheres. Slight deviations between the models for 
some quantities are probably caused by differences in spa- 
tial resolution and in the opacity binning procedures. These 
results give confidence in the reliability of comprehensive 
simulations as a tool for studying stellar atmospheres and 
surface convection. 
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